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ABSTRACT 


Active integral twist control for vibration reduction of helicopter rotors during forward 
flight is investigated. The twist deformation is obtained using embedded anisotropic 
piezocomposite actuators. An analytical framework is developed to examine integrally- 
twisted blades and their aeroelastic response during different flight conditions: fre- 
quency domain analysis for hover, and time domain analysis for forward flight. Both 
stem from the same three-dimensional electroelastic beam formulation with geometrical- 
exactness, and are coupled with a finite-state dynamic inflow aerodynamics model. A 
prototype Active Twist Rotor blade was designed with this framework using Active 
Fiber Composites as the actuator. The ATR prototype blade was successfully tested 
under non-rotating conditions. Hover testing was conducted to evaluate structural 
integrity and dynamic response. In both conditions, a very good correlation was ob- 
tained against the analysis. Finally, a four-bladed ATR system is built and tested to 
demonstrate its concept in forward flight. This experiment was conducted at NASA 
Langley Transonic Dynamics Tunnel and represents the first-of-a-kind Mach-scaled 
fully-active-twist rotor system to undergo forward flight test. In parallel, the impact 
upon the fixed- and rotating-system loads is estimated by the analysis. While discrep- 
ancies are found in the amplitude of the loads under actuation, the predicted trend 
of load variation with respect to its control phase correlates well. It was also shown, 
both experimentally and numerically, that the ATR blade design has the potential for 
hub vibratory load reduction of up to 90% using individual blade control actuation. 
Using the numerical framework, system identification is performed to estimate the 
harmonic transfer functions. The linear time-periodic system can be represented by 
a linear time-invariant system under the three modes of blade actuation: collective, 
longitudinal cyclic, and lateral cyclic. A vibration minimizing controller is designed 
based on this result, which implements classical disturbance rejection algorithm with 
some modifications. The controller is simulated numerically, and more than 90% of 
the 4P hub vibratory load is eliminated. 

By accomplishing the experimental and analytical steps described in this thesis, 
the present concept is found to be a viable candidate for future generation low- 
vibration helicopters. Also, the analytical framework is shown to be very appropriate 
for exploring active blade designs, aeroelastic behavior prediction, and as simulation 
tool for closed-loop controllers. 
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Chapter 1 


Introduction 


1.1 Background 

Rotorcraft has been a very important mode of aerial transportation due to its capa- 
bility of vertical take-off and landing, enabling many unique missions such as rescue 
operation at sea. However, it has also been under several serious constraints such 
as poor ride quality due to high levels of vibration [1] and noise, restricted flight 
envelope, low fatigue life of the structural components, and high operating cost. The 
primary source of those problems is the complex unsteady aerodynamic environment 
which is generated near the rotor blades mainly during forward flight [2]. An in- 
stantaneous asymmetry of the aerodynamic loads acting on the blades at different 
azimuth location is developed, and such asymmetry becomes more and more adverse 
as the forward flight speed increases. Therefore, the rotor system is the major com- 
ponent from which helicopter vibrations originate, and the resulting vibratory load 
becomes a dominant factor of reducing the life of fatigue-critical components and 
poor ride quality. These vibrations also limit the performance of the helicopters such 
as forward flight speeds, and tend to decrease payload due to the addition of extra 
vibration-alleviation devices. The coupling between the structural and mechanical 
components such as rotor, fuselage, engine, and transmission adds another degree of 
complexity to this problem. 

A typical aerodynamic environment of the helicopter main rotor during forward 
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vibrations. 


1.2 Helicopter Vibration Reduction 

There have been considerable efforts to reduce the vibration in helicopter [4, 5), 
and vibration alleviation methodologies employed by the helicopter designers may 
be categorized into the following three groups: 

1. Varying passive structural properties of the rotor system or fuselage by tuning 
its dynamic characteristics [6]; 

2. Employing passive or active vibration absorbing devices either at the rotating 
system or the fixed system [7]; 

3. Direct modification of the excitation forces, principally aerodynamic forces to 
reduce vibration. 

The first and second categories involve installation of vibration absorbing devices 
which produce counteracting inertial and damping forces. They are still used in most 
of the rotorcraft flying today although they also bring unavoidable penalties in terms 
of weight and tend to affect vibrations only at discrete points. Therefore, an effort to 
modify directly the excitation forces has been sought by the helicopter community, 
that is, to eliminate or reduce vibrations by modifying unsteady aerodynamic forces 
acting on the rotor blades. During the last two decades, this has been investigated 
and different implementations were attempted. Higher harmonic control (HHC) and 
individual blade control (IBC) are the typical examples of these efforts. Higher har- 
monic control is accomplished by manipulating a conventional swashplate to enable 
blade pitch control of a higher multiple frequency than an integer multiple of rotating 
frequency, i.e., (kb i l)/rev. Individual blade control installs a feathering actuator in 
each blade rather than modulating the swashplate, and allows for blade pitch control 
at arbitrary frequencies. Several outstanding results were obtained regarding vibra- 
tion reduction capability of these concepts, and they comprise of analytical studies 
searching for an optimal control scheme [8, 9, 10], wind tunnel tests with either small 
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1.3 Previous Work Related with Integral Twist 


Actuation 

1.3.1 Actuators Applicable for the Integral Concept 

An anisotropic actuator is required for the implementation of integral blade twist 
actuation with certain required characteristics. First, it must be flexible enough to 
be inserted in the curved shape of the blade assembly. Also, it should have its own 
structural integrity to withstand the pressure applied during blade fabrication and the 
external loads during the blade operation. It must have high levels of strain-inducing 
capability at an appropriately applied electric field. Anisotropy of the actuation is 
required so that tailoring in the blade design may be possible. Finally, cheap actuator 
is preferred considering the final blade cost. Candidates which are presently available 
include Active Fiber Composites (AFC) and Macro-Fiber Composites (MFC). 

The AFC is an anisotropic, conformable actuator, which can be integrated with 
the passive structure [34]. It was originally developed at MIT and now being com- 
mercialized by Continuum Control Corporation, Billerica, Massachusetts. The AFC 
actuator utilizes interdigitated electrode poling and piezoelectric fibers embedded in 
an epoxy matrix. This combination results in a high performance piezoelectric ac- 
tuator laminate with strength and conformability characteristics much greater than 
that of a conventional monolithic piezoceramic. Fig. 1-2 shows AFC packs being 
inserted in the blade assembly conducted as part of this thesis. However, some dis- 
advantages are also identified for this actuator: difficulty of processing and handling 
expensive piezoelectric fibers during actuator manufacturing and high actuator volt- 
age requirements. Basic material characterization and proof of concept of an integral 
twist-actuated rotor blade was investigated at MIT’s Active Materials and Structures 
Laboratory [35]. 

The MFC has been recently developed at NASA Langley based on the same idea as 
the AFC in using the piezoelectric fibers under interdigitated electrodes [36] . In this 
actuator, shown in Fig. 1-3, the piezoelectric fibers are manufactured by dicing from 
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[30], built and tested a 6-ft diameter two-bladed Froude-scaled rotor model with banks 
of piezoceramic crystal elements in ±45° embedded in the upper and lower surfaces 
of the test blade [31, 37]. Using dual-layer actuators, the active blade achieved 0.5° 
of maxim um experimental tip twist actuation, still lower than the 1° to 2° necessary 
for the possible vibration control applications. Bernhard and Chopra studied another 
twist concept that incorporates an active bending-torsion actuator beam within the 
blade spar [28], Tip twist angles of 0.15° to 0.5° (at 100 V) were achieved during 
Froude-scaled blade hover test. 

The most relevant work for this thesis, however, is the one conducted by Rodgers 
and Hagood [32] as part of a Boeing/MIT program sponsored by DARPA [38]. They 
manufactured and hover tested a 1 /6th Mach-scaled CH-47D blade in a two-bladed 
rotor where the integral twist actuation was obtained through the use of AFC . In order 
to design the blade structure and predict the actuation performance, a rudimentary 
single-cell active composite beam model [39] was used. Also, an intentional reduction 
by 50% on the baseline torsional stiffness was imposed and regarded to improve twist 
actuation. Hover testing on the MIT Hover Test Stand Facility demonstrated tip 
twist performance between 1° and 1.5° in the rotating environment. Boeing/MIT 
continues this work that eventually should lead to forward flight wind-tunnel tests 
and full-scale blade section manufacturing [38]. 

1.3.3 ATR Blade - Previous Work 

Another example of an integral blade twist concept has been studied by the author 
and his co-workers [33, 40] as part of a NASA/ Army/MIT Active Twist Rotor coop- 
erative agreement program. The structural design of the ATR prototype blade em- 
ploying embedded AFC actuators was conducted based on a newly developed analysis 
for active composite blade with integral anisotropic piezoelectric actuators [41]. The 
formulation is one of the first attempts for asymptotically-correct analysis of active 
multiple-cell beams presented in the literature. The approach is based on a two-step 
solution of the original three-dimensional electromechanical blade representation by 
means of an asymptotical approximation: a linear two-dimensional cross-sectional 
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twist actuation was experimentally observed of 1.1° peak-to-peak. 

Even though the new vibration reduction approach using the active materials 
technology showed promising results from the conceptual point of view [47 , 40] and 
successful preliminary hover testing with small-scaled models [32, 48], experimental 
forward flight tests were not performed until recently [49]. Moreover, an active aeroe- 
lastic environment to design, analyze, and simulate the behavior of integrally-twisted 
active rotor systems needed to be created. This should impact the ability to design 
active rotors and will support the design of their control law. 


1.4 Present Work 

This thesis concentrates on the study of vibration reduction using integral twist actu- 
ation on helicopter rotor blades. It includes the development of an analytical frame- 
work to identify dynamic characteristics of active rotor system in different flight 
conditions. The analytical model for forward flight combines a geometrically-exact 
theory for the dynamics of moving beams with active materials constitutive relations 
and a finite-state dynamic inflow theory for helicopter forward flight aerodynamics. 
The solution of this is conducted in time domain. Numerical results from the ana- 
lytical models are correlated with experimental data obtained from bench, hover and 
forward flight testing. For this purpose, two main steps are taken. First, a hover 
testing using the ATE prototype blade previously manufactured is conducted, and 
the numerical results from the analytical hover model are correlated with the ex- 
perimental data. Among other things, this step ensures the adequacy of the blade 
design and evaluates the twist actuation performance in the rotating condition. At 
the same time, the validity of the proposed analytical framework is demonstrated in 
both non-rotating and rotating conditions. Secondly, a four-active-bladed rotor sys- 
tem based on the prototype blade design is manufactured and tested in the open-loop 
control manner for the forward flight condition in the wind-tunnel. In parallel, an 
active aeroelastic analysis model is upgraded to deal with active twist rotor system 
during forward flight. Regarding the experimental work, an approach similar to the 
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helicopter blades 

2. Extend it to an aeroelastic analysis framework to simulate helicopter flight with 

active twist rotor blades and evaluate its response functions 

3. Correlate the analytical model with experiments conducted in the bench, hover 

and forward flight 

4. Assess helicopter vibration reduction capability using active twist rotor blades 

with an appropriate closed-loop control algorithm. 
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Chapter 2 


Analytical Framework 

2.1 Introduction 

For analyzing helicopter blades with embedded strain actuators, a framework is 
needed such that the effects of the active material embedded in the structure are 
kept throughout all the stages of the analysis. The framework should also contain an 
appropriate aeroelastic analysis component to predict the blade behavior under actu- 
ation during different operating conditions such as non-rotating, hover, and forward 
flight. For this purpose, a specific analytical framework for an active helicopter rotor 
system is proposed and its schematic is illustrated in Fig. 2-1. 

A base element from which the framework originates is the structural model of a 
general composite beam with embedded anisotropic actuators, and this corresponds 
to the dashed block at the upper part in Fig. 2-1, designated as “3-D electroelastic 
beam.” In this structural model, an asymptotical analysis takes the three-dimensional 
electromechanical problem of a rotor blade and reduces it into a set of two analyses: 
a linear analysis over the cross section and a nonlinear analysis of the resulting beam 
reference line. Such separation of the blade problem makes it convenient for aeroe- 
lastic analysis and consistently accounts for the active material effects. Using the 
structural model previously developed in [41, 45], the ATR prototype blade design 
was conducted, and its static twist actuation performance was evaluated [45]. Impor- 
tant elements included in the structural model are recapitulated in this chapter. 
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The previous structural model is extended and improved in this thesis to investi- 
gate the dynamic characteristics of the active rotor system in the bench, hover, and 
forward flight conditions. This corresponds to the lower dashed block in Fig. 2-1, 
called “Aeroelastic solution.” In this extension, the global beam analysis element 
from the previous structural model is combined with an appropriate aerodynamics 
model to compose an aeroelastic system. However, the system is still dependent upon 
the cross-sectional analysis regarding the beam cross-sectional properties. Solution of 
the resultant aeroelastic system is obtained in either frequency or time domain, and 
it includes the blade loads, hub vibratory loads, and blade motion. 

Specifically, for hover calculation, a mixed form of the geometrically-exact beam 
analysis model [44] is modified to account for integral actuation, and combined with 
finite-state dynamic inflow unsteady aerodynamics [50]. Based on the idea of small 
perturbation from a steady-state equilibrium position, frequency response functions 
of the active blade are determined using Laplace transform of the state-space repre- 
sentation. By performing hover analysis, more insight into the fidelity of the present 
modeling can be explored. Also, the behavior of the integral actuator under rotating 
condition can be assessed before examining the overall actuation authority of the 
ATR system. 

In forward flight analysis, the same geometrically-exact beam formulation is uti- 
lized but in displacement-based form. Also, the forward flight version of the same 
aerodynamics model is used, and their solution is performed in time domain. For im- 
plementation, an existing multi-body dynamics code [51] is modified for the needed 
active beam analysis. Time domain integration is selected since it is adequate for 
simulation of the blade response under open-loop actuation. This enables system 
identification for the sake of modem control, and ultimately closed- loop performance 
of ATR systems can be studied. 
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Figure 2-2: Two-cell thin- walled cross section beam 


Associated with this displacement field, the beam constitutive relation which relates 
beam generalized forces (axial force, twist, and two bending moments, respectively) 
with beam generalized strains (axial strain, twist curvature, and two bending cur- 
vatures) and corresponding generalized actuation forces is obtained in the following 
form: 
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where [tfy] is the stiffness matrix function of geometry and material distribution at 
the rotor cross section. 71 is the axial strain, is the elastic twist, and k 2 , are 
two bending curvatures. The actuation vector is a function of the geometry, material 
distribution, and applied electric field. Detailed expressions for the stiffness matrix 
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dependence is assumed linear here, with a linearization of the actuators properties 
conducted around the operating condition. In practice, some nonlinear behaviors of 
the ATR prototype blade with respect to applied voltage were observed during bench 

test [33]. 


2.3 Global Beam Analysis 

2.3.1 Mixed Form for Hover Analysis 

The nonlinear one-dimensional global analysis considering small strains and finite 
rotations is presented here as a direct expansion of the mixed variational intrinsic 
formulation of moving beams originally presented by Hodges [44], and implemented 
by Shang and Hodges [54]. The notation used here is based on matrix notation and 
is consistent with the original work of Hodges [44]. Some steps of the original work 
are repeated here to help understanding the modifications required in this extended 
active formulation. 

As shown in Fig. 2-4, a global frame named a is rotating with the rotor, with its 
axes labeled a 1? a 2 and a 3 . The undeformed reference frame of the blade is named b, 
with its axes labeled 6i, b 2 and b 3 , and the deformed reference frame named 5, with 
its axes labeled B\, B 2 and B 3 (though not shown in Fig. 2-4). Any arbitrary vector 
V represented by its components in one of the basis may be converted to another 
basis like 

V b = C^Va, V B = C Ba V a (2-4) 

where C ba is the transformation matrix from a to b, and C Ba is that from a to B. 
There are several ways to express these transformation matrices. C ba can be expressed 
in terms of direction cosines from the initial geometry of the rotor blade, while C Ba 
contains the unknown rotation variables. 

The variational formulation is derived from Hamilton’s principle which can be 

written as 

f‘[6(K -U)+W}dx l dt = SK (2.5) 

Ji i Jo 
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where u a is the displacement vector measured in the a frame, 6 is the rotation vector 
expressed in terms of Rodrigues parameters, &\ is the unit vector [1, 0, 0] , A is the 
3x3 identity matrix, v a and w a are the initial velocity and initial angular velocity 
of a generic point on the a frame. ( ) is a derivative with respect to time, and ( ) is 
a derivative with respect to the spanwise curvilinear coordinate. ( ) operator applied 
to a column vector is defined as ( ) mn = — e nm k{ )fc> with e nm * being the permutation 
tensor. 

To form a mixed formulation, Lagrange’s multipliers are used to enforce the sat- 
isfaction of the kinematical equations, Eq. (2.7). 

Manipulating the equations accordingly [54], one can obtain the a frame version of 
the variational formulation based on exact intrinsic equations for dynamics of moving 
beams as 

r SU a dt = 0 (2.8) 

Jt i 

where 


<5II a = J l {6u'?C T C ab F B +6uZ[{C T C ab P B y +ti a C T C ab P B ] 
+0?C T C ab M B ~ ^C T C o6 (eI + 7 )F B 
+8$l[{C T C ab H B y + u> a C T C ab H B + C T C ab V B P B ] 
-lF T a [C T C ab {e x + 7) - C*ei] -6F?u a 

-m T a { a + 1 + d -^)c ab K - m’ a e 

+SP T a (C T C ab V B -v a - u a u a ) - 5P T a u a 
A - J + ^CFC^SIb - w.) 

-JH T a 6 - Su^fa - 6iplm a }dxi 

- I (Sul Fa + W T a M a - IF T a Ua - SM T jI ( 2 . 9 ) 
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following transformation and interpolation are applied within each element [54] 


x = Xi + £A k, dx = A k d£, ( )' = ^ ) 


(2.13) 


8 XI a 

= <5ui(l — 

0 + 8v,i + i£,, 

U a = Ui 

H'a 

= Mi( 1 - 

0 + 8xpi +1 €, 

e = 6i 

8F a 

= 8Fi( 1- 

- £) + 8Fi+ i£, 

F b = Fi 

8M a 

= 8Mi( 1- 

— 4 ) + 8Mi + i£, 

M b = Mi 

6P a 

= 8Pi, 

P B = Pi 


8H a 

= Mi, 

H B = Hi 



where u„ 6„ F t , Mi, Pi and H l are constant vectors at each node i, and all 8 quantities 
are arbitrary. £ varies from 0 to 1. 

With these shape functions, the spatial integration in Eq. (2.12) can be performed 
explicitly to give 

£) {8uJf Ui + Ml hi + bF T i f Fi + 6M l f Mi + spj f Pi + 6Hj f Hi 

+8uJ +1 f ui +1 + Ml+iUi+i + 8Fj +1 f F . +1 +8M i+1 fM i+ i} 

= <5u^) + 1 jFjv+i + Mn+iMn+i — 8F n + iUn+i ~~ 8M 

—Suf F\ — dip^Mi 4 - SF^ U\ — 6Mi 6\ ( 2 - 14 ) 

where the ,. . ., f M{+1 are the element functions explicitly integrated from the 

formulation. 

In each element function, 7 and k should be replaced with an expression that is 
function of F B and Mb using the inverse of Eq. (2.11), along with the piezoelectric 
forcing vector F { B a) and M ( B ] . So do V B and Q B with an expression function of P B 
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Undeformed Stale 



Unwarped Cross Section 


Figure 2-5: Beam in the undeformed and deformed configurations 


during forward flight. The integral actuation forces and moments existing inside the 
blade structure are realized in the form of finite element loads to the passive beam in 
the modified time domain analysis. In what follows, an overview of the modifications 


to the formulation for the forward flight analysis is presented. 
The kinetic and strain energies of the beam are 


K 


-t 

2 Jo 


V B 

O.B 
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Pb 

H b 


> dx i 
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1 
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7 B 
kb 


r 


Fb 

M b 


dx i 


(2.17) 


The velocity-displacement and strain-displacement, relationships are expressed as 


’ Vs ' 


C baT C BbT u 

k J 

[ = 

C baT C BbT u> 


(2.18) 
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Substituting Eqs. (2.19) and (2.11) into (2.21), and integrating by parts yield the 
governing equations as follows 

(< c Bb c ba p B )' - (( o Bb c ba F B y = q + (c Bb c ba F^ ) y ( 2 . 22 ) 

(C Bb C ba H B y - uC Bb C ba P B - (C Bb C ba M B y + (u'~u ') T C Bb C ba F B = 

q - {C^C^Mfy + (u'Tu 'fc^C^F^ 

where q r = [_q T , q T J are the externally applied loads per unit span, measured in 
a. As described earlier, the effect of the actuation forcing vector is treated as an 
additional external load in the right-hand side of Eq. (2.22), while considering the 
tr ansf ormation between the inertial frame and the deformed configuration. 


2.4 Aerodynamic Analysis 

2.4.1 Hover Aerodynamics 

The external loads f a and m a along the B frame (Eq. 2.9) can be written as: 


Ib 


m B 


1 ) 

gPooCa < 


0 

(W B3 - f fii )W B3 - c -fWl > 
(ffil - W B z)Wb 2 - f Vb 3 + fgfii ^ 



a< 


1^32^1 - Vb 3 4- 
0 
0 


> 


(2.23) 


which is based on a thin airfoil theory [57], Here, Poo is the air mass density, c is the 
blade chord length, a is the lift curve slope, c do is the profile drag coefficient, Wb 2 
and W Bs are the components of the relative wind velocity in B frame, and fii is the 
component of the blade rotational speed along Bi direction. By the transformation 
between the frames ci and B . these forces and moments are converted to those in 
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(2.26) 


where 


BZ 


= /^; v /(2„ + i)(2t+i) 


m 


Hr 


n — 1 

E 

g=m,m+2,... 


//r 


29 + 1 


9 {t-q)(t + q+ 1) 


[A™]: diagonal matrix with K™ = 
[V n m ] : diagonal matrix with 


V m = < 

v n ' 


\/3 |a°| for (m,n) = (0,1) 


2\/3 | a? | otherwise 


where the absolute value is added to ensure the symmetry about the state of zero 
inflow, or about the zero thrust level. The right-hand side of Eq. (2.26) is regarded 
as a pressure integral. 


2.4.2 Forward Flight Aerodynamics 

The same finite-state dynamic inflow aerodynamics model presented in the previous 
section is also used for forward flight analysis. This aerodynamic theory was originally 
developed for both hover and forward flight conditions [56]. Moreover, the forward 
flight part of this model was already implemented in DYMORE. 

This model was constructed by applying the acceleration potential theory to the 
rotor aerodynamics problem with a skewed cylindrical wake. More specifically, the 
induced flow at the rotor disk was expanded in terms of modal functions. As a result, 
a three-dimensional, unsteady induced-flow aerodynamics model with finite number 
of states was derived in time domain. This model falls on an intermediate level 
of wake representation between the simplest momentum and the most complicated 
free wake methodologies. It does not require a severe computational effort, which 
is usually the case in those that involve the vortex filament theory. Therefore, this 
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can be expressed in the form 


'X ' 


' X' 


’X(t)' 


> = i 


> + < 






. . 


(2.30) 


where ( ) denotes steady-state solution and ( ) denotes the small perturbation about 
it. 

For the steady-state solution, one has to solve a set of algebraic nonlinear equations 
originated from Eq. (2.27) and Eq. (2.29): 


F s (X,0,V)-F L (X,Y,0) = 0 
-Fp(X,Y) + Fi(Y,0)=0 


The Jacobian matrix of the above set of nonlinear equations can be obtained 
analytically and, even with the modifications caused by the active material embedded 
in the structure, it is found to be very sparse. Note that the presence of the actuation 
in the blade changes the original terms of the Jacobian in a similar manner it does 
in Eq. (2.14). The steady-state solution can be found very efficiently using Newton- 


Raphson method. 

In order to investigate the dynamic response of the blade with respect to voltage 
applied to the embedded anisotropic strain actuator, a state-space representation is 
required once the steady-state solution is obtained. Perturbing Eqs. (2.27) and (2.29) 


using Eq. (2.30) about the calculated steady state yields 


r flEs. _ Ml 9Fl 

ax ax by 

aF P dfj _ aFp 

— ax dY BY 



’ X' 


4 


■ + 

X=X 




S t 



QFp _ dF*. 

ax ax 


BF, 


_ 

’ X' 

< 

> 

x=x 



[y J 


+ 




Y=Y 


[Ms. ol 


'O' 

av u 

{V} = < 

> 

0 0. 

■>< 

II 

X 



II 

.0, 


(2.32) 


57 


and all other zeros. Eq. (2.35) becomes the second part of a state-space representation 
with a generalized form as 

y = CX + T>V (2.36) 

where y is the output vector corresponding to the sensors embedded along the blade. 

Frequency response function of the blade can be calculated using Laplace trans- 
form of the simultaneous equations which are composed of Eqs. (2.33) and (2.36) 

HQ = C(Es- A) _1 B + D (2.37) 

V{*) 

Note that the coefficient matrix E is usually singular due to the mixed formulation 
of the beam model. 

2.5.2 Time Domain Solution for Forward Flight Analysis 

While the hover analysis presented in the previous section seeks frequency-domain 
quantities of the blade response function, the forward flight analysis performs a di- 
rect time integration of the blade response due to an integral actuation. This is 
needed since system identification and open- and closed-loop simulations, all in time, 
must be conducted. DYMORE, the original passive blade dynamics model, adopts 
a time-discontinuous integration scheme with energy decaying characteristics in or- 
der to avoid high frequency numerical oscillation [51, 58]. Such a high frequency 
oscillation usually occurs during a finite element time integration of a complex multi- 
body dynamic system. Details of the energy decaying time integration of the beam 
formulation are found in [51, 58], and briefly summarized in Appendix B. 

DYMORE is also capable of adjusting its time step size automatically to maintain 
stability and accuracy of the integration scheme. Another advantage of adopting 
multi-body formulation here is that the total shear force and moment exerted by the 
rotor system can be easily extracted. By adding and monitoring a rigid body element 
which represents a rotor shaft, the degree of vibratory load variation of the entire rotor 
system can be directly evaluated. Finally, the control sensitivity functions due to high 
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Chapter 3 


Experimental Setup 


3.1 Overview 

In the previous chapter, an analytical framework is proposed and established to pre- 
dict aeroelastic behavior of an active rotor system and numerically evaluate its effec- 
tiveness in vibration reduction. At the same time, an experimental effort was pursued 
to substantiate the present integral blade actuation concept through a small-scale 
wind-tunnel model. Results from these experiments are also utilized for correlation 
with the predictions from the proposed analytical framework. Wind-tunnel tests were 
conducted at NASA Langley’s Transonic Dynamics Tunnel as part of the collabora- 
tion between the U.S. Army Research Laboratory, at NASA Langley Research Center, 
and MIT. 

The ATR prototype blade was previously designed and successfully manufactured. 
Preliminary bench testing was conducted to confirm its basic structural characteris- 
tics. Details of the relevant work are found in [45], and summarized in this chapter. 
The prototype blade is used for hover test with three other dummy blades to compose 
the four-bladed fully-articulated rotor. Blade response under rotating condition is in- 
vestigated in the hover test. A minor modification is added to the prototype blade 
design, and four active blades are fabricated based on it. Using four-active-bladed 
rotor system, forward flight test is performed in an open-loop control manner. The 
results from the hover, and forward flight tests will be correlated with those predicted 
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(Unit:inch) 


Figure 3-1: Planform and cross-section of the ATR prototype blade (Dimensions are 
in inches.) 



Rotor Speed (RPM) 

Figure 3-2: Fan plot of the ATR prototype blade from the proposed analysis 
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prototype blade was originally expected to achieve static twist actuation amplitudes 
of between 2.0° to 2.5° and hovering flight dynamic twist actuation amplitudes of 
2.0° to 4.0° (based on CAMRAD II and PETRA simulations [40]) at the extended 
cycle of maximum applied voltage of 4,000 Vp P /l,200 Vpc- Structural integrity of 
the new blade design was evaluated based on the worst loading conditions, which are 
expected to occur within the rotor system operating envelope. In this design, forward 
flight with the maximum speed is selected as the design loading condition. Then, the 
largest magnitudes of the aerodynamic loads are extracted and combined with the 
centrifugal loads in order to give the worst loading values. A safety factor of 1.5 was 
used. 

3.2.2 ATR Test Blade with Modification 

Even though experimental structural characteristics and twist performance of the 
prototype blade compared well with design goals [45] , a concern was raised regarding 
its structural integrity, especially affected by the fatigue loadings. There was not 
enough experimental evidence that the prototype blade had the fatigue life according 
to the criterion employed by NASA Langley Research Center for wind-tunnel testing 
models [59]. Also, an empirical formula adopted by the contractor who planned 
to build the ATR test blades for forward flight testing indicated that the design 
should be improved in fatigue [60]. A modification was applied to the design of the 
prototype blade to compensate this shortcoming. Different lay-up configurations were 
suggested to increase the structural integrity within the range that its characteristics 
is not greatly changed from that of the prototype blade. As a result, only one ply 
of E-Glass fabric prepreg in 0/90° was added to the front spar assembly in order 
to further withstand centrifugal loading. Using the active cross-sectional analysis 
previously described, basic structural characteristics were computed for the updated 
configuration, and listed in Table 3.3. The material properties of the passive prepregs 
and the AFC plies used in the test blades manufacturing are slightly changed from 
those in the prototype blade, and are summarized in Appendix D. 
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Figure 3-3: ATR prototype blade 


3.4 Aeroelastic Tests 

3.4.1 Wind Tunnel 

The Langley Transonic Dynamics Tunnel (TDT), whose schematic is shown in Fig. 3- 
4, is a continuous-flow pressure tunnel capable of speeds up to Mach 1.2 at stagnation 
pressures up to 1 atm. The TDT has a 16- ft square slotted test section that has 
cropped corners and a cross-sectional area of 248 ft 2 . Either air or R-134a, a heavy gas, 
may be used as a testing medium. The TDT is particularly adequate for rotorcraft 
aeroelastic testing due to several advantages associated with the heavy gas. At first, 
the high density of the testing medium allows model rotor components to be heavier, 
and this satisfies the structural design requirements easily while maintaining dynamic 
scaling. Second, the low speed of sound in R-134a (approximately 170 m/sec) allows 
lower rotor rotational speeds to match full-scale hover tip Mach number. Finally, the 
high-density environment increases Reynolds number throughout the testing envelope, 
which enables more accurate modeling of the full-scale aerodynamic environment of 
the rotor system. Both hover and forward-flight tests of the ATR system are primarily 
conducted in the heavy gas testing medium at a constant density of 2.432 kg/m 3 . 

3.4.2 Test Apparatus 

The Aeroelastic Rotor Experimental System (ARES) helicopter testbed, whose sche- 
matic drawing is illustrated in Fig. 3-5, is used for both hover and forward-flight 



1.35 


- 8.37 



Figure 3-5: Schematic of the Aeroelastic Rotor Experimental System (ARES) heli- 
copter testbed (All dimensions are in ft.) [48] 

the fixed-system forces and moments generated by the rotor model. A streamlined 
fuselage shape encloses the rotor controls and drive system. However, the fuselage 
is isolated from the rotor system such that fuselage forces and moments do not con- 
tribute to the loads measured by the balance. 

Fig. 3-6 shows the ATR blades mounted on the ARES helicopter testbed in the 
TDT. For this configuration a four-bladed articulated hub with coincident flap and lag 
hinges is used on the ARES. The feathering bearing for the hub is located outboard 
of the flap and lag hinges, and trailing pitch links are used. The hub is configured 
such that pitch-flap coupling of 0.5 (flap up, leading-edge down) is obtained and the 
lag-pitch coupling is minimized. During the hover testing, the test section floor, and 
the ARES testbed, is lowered approximately 3 ft to allow the rotor wake to vent into 
the surrounding plenum volume, thus reducing recirculation effects. (In Fig. 3-6, the 
test section floor is shown in its normal, raised position.). 
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Table 3.4: Hover test conditions for the ATR prototype blade 


Testing 

Density 

Rotor speed 

Collective 

Voltage 

medium 

(kg/m 3 ) 

(rpm) 

pitch (deg) 

amplitude (V) 

Air 

1.225 

400 

0 1 

100 



400 

0 

500 



400 

0 

750 



400 

0 

1000 

Air 

1.225 

688 

0 

500 



688 

0 

1000 



688 

4 

500 



688 

4 

1000 



688 

8 

500 



688 

8 

1000 



688 

12 

1000 

Heavy 

2.432 

688 

0 

500 

gas 


688 

0 

1000 


688 

4 

500 



688 

4 

1000 



688 

8 

500 



688 

8 

1000 

Heavy 

1.546 

688 

8 

1000 

gas 

1.984 

688 

8 

1000 

2.432 

688 

8 

1000 


2.432 

619 

8 

1000 


to 100 Hz, in 5 Hz increments, at amplitudes of up to 1000 V were applied to the 
ATR prototype blade. Data from the blade strain gauge bridges, the ARES testbed, 
and the high-voltage amplifier channels were recorded at a rate of 3,000 samples-per- 
second by the computer control system for 5-second durations. The signals acquired 
through the channels during the test are listed in Table 3.5. Even though different 
loads were measured in the test, correlation with the analytical framework is limited 
to the blade torsion moments in this thesis. The total rotor loads were measured at the 
fixed-system balance, and since the dynamics of the fuselage model (housing motors, 
slip rings for data and power, etc.) is not completely available, these quantities are 

not included in the present hover analysis model. 

The acquired signal data in time were processed to obtain the transmissibility of 
the system with respect to the sinusoidal actuation of the active blade in frequency 
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Table 3.6: Forward flight test conditions for the ATR system 
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sine-sweep signal was applied by the high-voltage amplifiers. In case of sine-dwell 
signals, only 3P, 4P, 5P frequency components were considered since (6+1, 6, 6-1) 
frequency components are to influence significantly 6-bladed rotor system. Available 
blade control modes include collective twist, differential twist, and an Individual Blade 
Control (IBC) mode where each blade actuates according to a prescribed schedule 
with respect to its position in the azimuthal location. In collective twist mode, all 
the blades are under the same synchronous twist actuation signals, while those of an 
opposite sign are transferred to the blade at opposite azimuthal location (e.g., Blade 
No. 1 and 3) in differential mode. For IBC actuation mode, the actuation on each 
blade behaves in the same phase at a specific azimuthal location. 

Also, a sweeping algorithm over control phase angle was considered within the 
IBC scheme. It is worth noting that although control phase is indicated by the rotor 
azim uth (i.e. , 0° control phase is coincident with 0° azimuth), control phase is not 
equivalent to rotor azimuth. For example, a 3P twist actuation with control phase of 
180° would impose the maximum twist control at a rotor azimuth of 60°. The second 
and third cycles would achieve maximum twist control at 180° and 300 , respectively. 
This will be revisited in the IBC signal generation for the forward flight analysis in 
Section 5.3.2. 

The sine-sweep signal was also used for experimental system identification pur- 
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Chapter 4 


Characteristics of the ATR Blade 
on the Bench and in Hover 


4.1 Overview 

Using the ATR prototype blade, two major tests are executed regarding its dynamic 
response induced by a time-varying electric field applied to the embedded AFC actua- 
tors. Bench top non-rotating actuation testing is conducted first. Then, the prototype 
blade is used to build an active rotor system with the other dummy blades to be tested 
in hover condition. Details of the hover experiments plan are already introduced in 
Chapter 3. All these test data are compared with the results from the analytical 
framework presented before, especially the model which obtains the solution of the 
aeroelastic system in frequency domain. 

4.2 Basic Bench Testing 

Different characterization tests were performed on the ATR prototype blade at bench 
top condition in order to validate the design and manufacturing procedures, and to 
verify the performance of the prototype article. Also numerical results obtained from 
the proposed analytical framework were compared with the test data at each stage 
of the experiments since the model was developed and expanded to enable variety 
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Table 4.1: Peak-to-peak tip twist actuation of the ATR prototype blade (2,000 V pp /0 
V DC , 1 Hz) 



Present 

Experiment 

Experiment 


Model 

(LDS) 

(PMI) 

Spar only 

1.4° 

1.1° 

- 

Spar+fairing 

1.2° 

1.0° 

1.1°* 


* only 18 active AFC packs. 


experimental result from the PMI test is about 15% higher than the original LDS 
measurements, and the difference from the analytical model is within 12%. This dis- 
crepancy was expected and arises from variation in AFC material properties (between 
packs) and uncertainty in the material properties used in the analysis. 


4.3 Non- Rotating Frequency Response 

The non-rotating dynamic characteristics of the prototype blade can be evaluated 
from the frequency response of an applied sinusoidal excitation to the AFC actuators. 
Figs. 4-1 - 4-4 show the results of the laser displacement sensor (for tip twist angle), 
the blade strain gauges readings, and the predictions of the frequency domain analysis 
at several blade stations as function of the AFC actuation. 

As one can see, the first torsional mode is clearly identified at approximately 85 
Hz, and this result matches well with model prediction. The model neither includes 
structural nor stationary aerodynamic damping, resulting in infinite peaks at reso- 
nance. This already indicates that some structural damping should be added to the 
model. Once the aerodynamics is included in the problem, its damping will bring that 
to a finite amplitude. While the strain results could be obtained at high excitation 
voltages, the dynamic tip twist was measured at low voltages due to the limitation 
on the range of the laser sensors. At 400 V pp , the peak-to-peak tip twist response 
of the blade is approximately 3.5°. Such an increased dynamic response around the 
first torsional natural frequency is expected to affect the twist response over the fre- 
quency range of interest when the blade is rotating. It makes the frequency response 
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Figure 4-3: Equivalent torsional moment at 49% blade radius of the ATR prototype 
blade on the bench 



Figure 4-4: Equivalent torsional moment at 75% blade radius of the ATR prototype 
blade on the bench 
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Figure 4-5: Equivalent torsional moment at 31% blade radius of the ATR prototype 
blade in hover (688 rpm, heavy gas, medium density = 2.432 kg/m 3 , 2,000 V pp actu- 
ation, e o = 0°, 4°, 8°). 


varying from 40% to 25%, respectively. This constant offset may indicate the effects 
of the local three-dimensional deformation field induced by the presence of collocated 
actuators on the strain gauge bridge, and its effects on strain gauge calibration not 
taken into account in these results. The phase component of the predicted equivalent 
torsional moment correlates very well with the experimental results, with errors less 
than 8% when approaching 10P, which is associated with initial saturation of the 
power amplifiers during tests. 


4.4.2 Medium density sensitivity 

When changing the testing medium density, the resulting frequency response func- 
tions are shown in Figs. 4-6 and 4-7 for equivalent torsional moment obtained at 31% 
and 51% spanwise locations, respectively. As one can see, the medium density vari- 
ation does not influence the actuation authority, except at the torsional resonance 
frequency due to the change in the aerodynamic damping with density. It is also 
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Figure 4-7: Equivalent torsional moment at 51% blade radius of the ATR prototype 
blade in hover (688 rpm, 2,000 V pp actuation, 0 O = 8°, medium density = 1.546 
kg/m 3 , 1.984 kg/m 3 , 2.432 kg/m 3 ). 


Model Root/Blade Channels/S:T1 BS17.0 




Figure 4-8: Equivalent torsional moment at 31% blade radius of the ATR prototype 
blade in hover (2,000 V pp actuation, 9 a = 8°, medium density = 2.432 kg/m 3 ). 
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Tip Twist Response 




Figure 4-10: Blade tip twist amplitude predicted by the proposed analytical frame- 
work (2,000 V pp actuation, 680 rpm, 9 0 = 8°, medium density = 2.432 kg/m 3 ) 

the experimental setup to avoid ground and air resonance of the rotor system. The 
coupled pitch-flap-lag motion may bring some of those effects to influence the results 
above. Forward flight part of the proposed analytical framework has a capabilty of 
modeling lead-lag damper, therefore better correlation is expected with regard to this 
matter. 

Since no specific sensor for tip twist measurements were included in the prototype 
blade, the blade tip twist can only be estimated based on the analytical framework. 
Fig. 4-10 presents such results. As it can be seen, between 3P and 5P, the blade tip 
twist amplitude varies between 1.0° and 1.3°, which, according to previous CAMRAD 
II simulations [40], should be enough to provide 60% to 80% reduction on 4P hub 
shear vibratory loads. Forward flight test and analysis in Chapter 5 will be addressing 
this issue in detail. 
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Chapter 5 


Dynamic Characteristics of the 
ATR System in Forward Flight 

5.1 Overview 

The forward flight regime is of great interest for the vibration reduction problem. 
The blade twist control is suggested to alter the undesirable unsteady aerodynamic 
environment which develops in that flight regime. As mentioned previously, wind- 
tunnel testing is conducted on the active rotor system with the ATR test blades. At 
the same time, confidence on the established analysis model for forward flight is to 
be obtained through its correlation with experimental data. 

Initially, the bench top static actuation testing is revisited in this chapter. This is 
to exemplify basic validation of the active time domain analysis. Then, the potential 
impact upon the fixed- and rotating-system loads by the integral blade actuation 
during forward flight is examined. Both experimental and analytical efforts focus on 
an open-loop control and their correlation. By accomplishing this, the present forward 
flight model can be taken for further analytical tasks related with system identification 
and closed-loop controller design, which will be introduced in Chapters 6 and 7. 
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Time (sec) 

Figure 5-1: Time history of tip twist angle of the ATR blade at bench by sine-sweep 
actuation signal 



Excitation Frequency (Hz) 


Figure 5-2: Tip twist response of the ATR blade at bench 



Active 



Figure 5-4: Detailed multi-body representation of 4-active-bladed ATR system 


to give the hub vibratory loads. Since the ATR system is fully-articulated, three 
revolute joints are consecutively located between the root retention and the active 
blade to represent flap, lead-lag, and feathering hinges. As shown in Fig. 5-4, the 
flapping and lead-lag hinges are coincident. Among the three joints, a prescribed 
collective and cyclic pitch control commands are applied at the feathering hinge, and 
their numerical values are based on those used in the wind-tunnel experiment. These 
are summarized in Table 5.1. Finally, active beams are attached to represent the ATR 
blades, and they are discretized during the analysis with at least four beam elements 
per blade, each with the 3rd-order interpolation polynomials. Therefore, there are 
approximately 900 degrees of freedom to be solved at each time step, including the 


Table 5.1: Trim control inputs for the forward flight test conditions 


Advance 

Rotor shaft 

Collective 

Longitudinal 

Lateral 

ratio 

inclination angle 

Pitch 

cyclic pitch 

cyclic pitch 


(deg) 

(deg) 

(deg) 

(deg) 

0.14 

-1.0 

7.5 

-3.5 

-3.1 

0.17 

-1.0 

7.0 

-3.6 

-3.0 

0.20 

-1.0 

7.1 

-3.9 

-3.1 

0.233 

-2.0 

7.5 

-4.3 

-3.3 

0.267 

-2.0 

7.8 

-4.8 

-3.4 

0.30 

-4.0 

8.0 

-5.0 

-3.3 

0.333 

-6.0 

10.7 

-6.2 

-4.2 

0.367 

-6.0 

11.2 

-6.7 

-4.4 
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HV signal history 



Figure 5-5: Example of high-voltage input generated for an IBC-mode 3P actuation 
with 12 divisions of control phase angle 


those under actuation. Then, for each 0.5-s period of actuation, each with different 
control phase angle, and another 0.5-s period of no actuation is applied between them. 
These are applied one after the other as shown in Fig. 5-5. 

In each 0.5-second period of actuation cycle, there exists a sine-dwell signal corre- 
sponding to a frequency 3P with different phase angles for each blade. 3P sine-dwell 
signals generated in this fashion are clearly seen in the magnified plot at the right 
side in Fig. 5-5. By applying this control phase algorithm, the blades exhibit the 
maximum amplitudes of the sinusoidal electric field at certain azimuthal locations as 
exemplified in Fig. 5-6. The maximum amplitude occurrence during the first and sec- 
ond actuation periods shows an azimuthal difference of 10° corresponding to 120°/12 
divisions, although it is designated as a phase difference of 30° corresponding to 
360°/ 12 divisions in terms of control phase. 

5.3.3 Results of the Model without Pitch Link 

As a result of the simulation, a time history of the quantities of interest, for ex- 
ample, the forces and moments exerted at each blade and root retention, blade tip 
displacements, flapping and lead-lag motions at the articulated hinges, and aerody- 
namic forces generated at the blades are recorded. Among the flight conditions tested 
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Figure 5-7: Simulated time history of hub vertical shear forces when the 3P actuation 
is applied as described in Fig. 5-5 


for 3P, 4P, 5P actuation applied during the same steady-state trim condition. In the 
figures, the lines are simple interpolation of the solution points which are obtained 
from the analysis at discretely increased control phase, while the experimental data 
are still displayed with the discrete symbols. 

The variation of the 4P hub vibratory load components are calculated with respect 
to the variation of the control phase angles in the vertical, forward, and sideward 
directions, respectively. The hub sideward component is not included in Fig. 5-8 
since its magnitude is too high compared to the other force components, indicating a 
problem with the measuring device. 

The load predicted from the analysis shows significant discrepancy in amplitude 
from the experimental results, although their variation trends in terms of control 
phase are in good agreement. As one can see from Fig. 5-8, 3P frequency sine-dwell 
actuation appears to be the most effective in reducing the hub shear vibratory loads 
in both cases of vertical and forward components. More specifically, 3P actuation 
is most effective in hub vertical shear load reduction, resulting in 95% reduction at 
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210° control phase. It also reduces the hub forward shear loads by 80% at about 180° 
control phase. Such a hub shear vibratory load reduction performance numerically 
predicted here shows similar trend as it was observed in the experiment. Further 
discussion on the comparison between the current analytical results and the experi- 
mental data will be given at the end of this chapter. Since the dynamics of the test 
apparatus (ARES) used for the wind-tunnel test is not included in the model, this 
may be responsible for the discrepancies. Upgraded input model for the same ATR 
system including the pitch link and all the linkage components in the swashplate is 
attempted for better correlation, and will be described in a later section. Also, experi- 
mental characterization of the ARES testbed used in the wind-tunnel test is expected 
to be conducted in the future for the precise modeling of these extra components. 

Rotating-system loads in the low-speed level flight case 

Quantities in the rotating frame, for example, the flap and chordwise bending mo- 
ments, and torsional moments are calculated in the reference blade and can also be 
correlated with the experimental results. The span location where these quantities 
are calculated is selected to match those of the strain gauges embedded in the test 
blade. While the fixed-system quantities, such as hub shear vibratory loads, were 
investigated only in 4P frequency components, those in the rotating frame are ex- 
tracted and examined in their 3P, 4P, and 5P frequency components. Figs. 5-9 - 5-10 
show blade loads at those frequencies. 

Again, the lines are simple interpolation of the solution points from the analysis, 
and symbols represent the experiments. The 3P frequency components of the flap 
bending moment at 28.7% span location are extracted and shown in Fig. 5-9 (a). The 
results are for the condition of /x = 0.140, as = —1°, Cr = 0.0066, and 1,000 V twist 
actuation at 3P, 4P, 5P with respect to control phase. 4P frequency components for 
the same flap bending moment are presented in Fig. 5-9 (b), and 5P components are 
in Fig. 5-9 (c). As well as in the fixed-system quantities, the present model captures 
the trend of variation in the rotating frame values as they were observed in the exper- 
iments. However, discrepancies can be observed in the amplitude by approximately 
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Figure 5-10: Variation of torsional moment at 33.6% span location for - 0.140, 
as — -1°, Ct = 0.0066, and 1,000 V twist actuation at 3P, 4P, 5P with respect to 

control phase 
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Figure 5-11: Polar plot of 4P hub shear vibratory loads for /z = 0.140, a s = -1 
C T = 0.0066, and 1,000 V twist actuation at 3P, compared with the baseline value 


01 



Control Phase (deg) 


(a) Vertical 



Control Phase (deg) 

(b) Forward 

Figure 5-12: Variation of 4P hub shear vibratory loads for p = 0.333, a s = -6°, 
C T = 0.0066, and 1,000 V twist actuation at 3P, 4P, 5P with respect to control phase 
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Figure 5-14: Variation of torsional moment at 33.6% span location for /x — 0.333, 
o; s = -6°, C T = 0.0066, and 1,000 V twist actuation at 3P, 4P, 5P with respect to 

control phase 
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Figure 5-15: Upgraded input model of ATR system including pitch link and swash- 
plate components 

flexibility of the whole control system. According to NASA Langley, the rigid feather- 
ing mode of the ATR system is observed at frequency between 15P and 16P. Based on 
the pitch inertia of the ATR test blade and pitch link geometry, an axial stiffness of 
the pitch link is estimated as also included in Table 5.2. By including the pitch link 
mechanism in the model, the pitch/flap coupling existing in the ARES apparatus, 
which amounts to 0.5 in flap up/leading edge down fashion, is automatically imple- 
mented. Due to the dynamic interaction among the control linkage, the feathering 
motion at the blade root in the rotating condition shows discrepancy from those in 


Table 5.2: Geometry and material property of the upgraded ATR system model 


Parameter 

Value 

6flap, lag 
^pitch 

XPL 

ypL 

Xsp 

ysp 

EA pitch link 

0.0762 m 
0.1143 m 
0.03556 m 
0.05715 m 
0.027686 m 
0.04445 m 
3.822 10 3 N 
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Figure 5-16: Variation of 4P hub shear vibratory loads for /x = 0.140, as = -1°, 
C T = 0.0066, and 1,000 V twist actuation at 3P, 4P, 5P with respect to control 
phase: experiment, analysis without pitch link, and analysis with pitch link 
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Figure 5-18: Variation of torsional moment at 33.6% span location for [i = 0.140, 
as = —1°, Ct = 0.0066, and 1,000 V twist actuation at 3P, 4P, 5P with respect to 
control phase: experiment, analysis without pitch link, and analysis with pitch link 
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5.4 Correlation of Forward Flight Analyses with 


Experiments 

Both of fixed- and rotating-system loads predicted from the forward flight analysis 
exhibit significant discrepancy in amplitude from the experimental results, although 
their variation trend in terms of control phase is in good agreement. It should be 
further noted that the baseline amplitude in most of the cases are significantly un- 
derpredicted by the numerical analysis. 

A study to improve the correlation of the baseline load amplitude is performed. 
For the analysis model without pitch link, blade structural modeling accuracy was 
studied with respect to the stiffness matrix and chordwise c.g. location, regarded 
as factors which affect the dynamics of the rotor system. According to the cross- 
sectional analysis results shown in Table 3.3, the ATR test blade has its chordwise 
c.g. at 18%, although the blade manufacturers suggested that the blade had required 
25% c.g. location. Also, the analysis indicates that the shear center is located at 
30% chord instead of the 25% chord that coincides with the blade reference line. This 
effect can be included in the analysis by simply using the fully-populated stiffness 
matrix. Therefore, different models of the ATR dynamics were constructed and their 
prediction of the baseline loads were obtained for the low-speed level flight condition, 
as shown in Table 5.3. Changes in the load prediction resulted from the variation 
in the structural modeling. However, these changes fail to significantly improve the 
correlation with the experiments. Finally, the model with pitch link did exhibit an 
improvement even with a crude structural modeling for the blade; however, it is 
a g ain a slight improvement. Therefore, an analysis model with a diagonal stiffness 
representation, 30% shear center and 25% c.g. locations, has been used for all' the 
analytical results presented in this chapter, and is also used for further analysis in 
the following chapters. 

Also, the experimental 4P hub sideward baseline force (See Table 5.3) presented 
unreasonably high magnitude when compared to the other two hub force components. 
This suggests that further characterization of the whole experimental apparatus is 
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Chapter 6 


System Identification of the ATR 
System in Forward Flight 

6.1 Overview 

The forward flight analysis model established in this thesis exhibited sufficient details 
of a typical helicopter and its rotor blade dynamic behavior, although it showed 
inaccuracy in amplitude of the predicted loads. Therefore, it is selected to be used in 
this chapter for system identification of the ATR system in forward flight. 

During forward flight, the helicopter rotor blade exhibits an aerodynamic envi- 
ronment which varies itself with a period corresponding to the rotor revolution. This 
situation is illustrated in Fig. 1-1. This signifies that the helicopter rotor system 
is basically a linear time-periodic (LTP) system during forward flight. Therefore, 
a methodology considering this periodicity is required for its characterization. In 
this thesis, a method is adopted which results in multi-component harmonic trans- 
fer functions [62]. The theoretical background of the adopted methodology is briefly 
summarized in Appendix E with its implementation schemes that include additional 
assumptions imposed on the transfer functions. The sine-sweep input signals cre- 
ated for computation of such harmonic transfer functions are described in detail in 
this chapter. The results of the system identification is presented for each different 
mode of blade actuation. Finally, certain characteristics of the present ATR system 
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equation., Eq. (6.2), will give the phase angle of the chirp <p c &S 

«0=(/ot + ^* 2 )2^ (6.3) 

where the phase is in radians. The sine-sweep signal is generated using the unwrapped 
azimuthal location quantities. The relation between rotor rotational speed ft (rpm) 
and the azimuthal location xp is 

(6.4) 

where time t is in seconds. A pseudo-time t is introduced based on the previously 
established relations, and uniformly distributed phase of N chirps over 360° can be 
produced by considering the number of the chirps that have already been generated, 
n c , and shifting the time vector accordingly, so that 

* = [^-{-360 + mod(360^,360)}]J^ (6.5) 

where “mod” is the modulo function that returns the remainder obtained from the 
division of two arguments. The constructed chirp signal vector U c with amplitude 
A c for the rotor system can be represented as 

U c = A c sin [27 t</> c (tjj V (6-6) 

where V is a vector of length b, and collective, cyclic, and differential mode of actu- 
ation among the blades can be achieved by adjusting the elements of V. In case 
of four-bladed rotor system, V for collective mode will be [ 1 1 1 1 } T , while 

j 1 ^ —l} T for differential mode. A block diagram which generates the sine- 

sweep input signals described so far is illustrated in Fig. 6-1, and it is a part of the 
actual Simulink program used to implement chirp actuation during the forward flight 
testing. 

Using the algorithm established, a collective mode chirp input signal is constructed 
with the amplitude of 1,000 V and nine phase angle divisions over 360° for the present 
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Figure 6-2: Collective mode sine-sweep signals generated with 9 divisions of phase 
angles 

four-bladed ATR system identification. According to the definition of collective mode, 
all four blades have the same synchronous input signals for actuation, therefore one 
signal corresponding to Blade No. 1 may become a representative of all the input 
signals. In Fig. 6-2, representative signal generated for each different phase angle 
is overlapped with one another along the azimuthal location. It can be clearly seen 
that the initiation phase angle for each signal is uniformly separated, therefore all the 
signals conduct one complete sweep of 360° azimuthal location. 

6.3 Results of the ATR System Identification 

6.3.1 Collective Mode of Actuation 

Using the constructed sine-sweep input signal, fixed- and rotating-system response 
of the ATR system at the low-speed level flight condition which was considered in 
Chapter 5 are calculated by the time domain analysis. At first, a series of collective 
mode actuation signal, as shown in Fig. 6-2, is applied, and its response is examined 


119 








Transfer Function G_ 

O 




Figure 6-6: Harmonic transfer function G 0 of the hub vertical shear loads during the 
collective mode actuation 

in Fig. 6-6. Furthermore, the property of G 0 being able to approximate the LTP 
system can be verified by comparing with the polar plot of the hub vertical shear 
loads obtained in Section 5.3.3. An amplitude of 47 N and phase delay of 317° can 
be read from the Bode diagram of Go in Fig. 6-6 at 4P excitation frequency. This 
information exactly matches those included in the polar plot of the hub vertical shear 
loads under the sine-dwell actuation at frequency 4P, which is illustrated in Fig. 6-7. 
In fact, IBC-mode 4P sine-dwell actuation considered in Chapter 5 generates exactly 
the same actuation signal as the collective mode considered here. The maximum 
amplitudes during IBC-mode 4P sine-dwell actuation occur at 90 , 180 , 270 azimuth 
after blade No. 1, respectively, which eventually places them at Blade No. 2, No. 
3, No. 4, respectively. This results in equivalent signal used in the collective mode 
of actuation. Therefore, both of open-loop control simulation using IBC-mode sine- 
dwell actuation in Chapter 5 and system identification using sine-sweep actuation 
here describe the system behavior in a consistent manner. 
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Figure 6-8: Matrix of G 0 estimated for three components of 4P hub shear vibratory 
loads versus three modes of blade actuation 


still valid for these modes of actuation. All the identification results including the 
collective mode previously obtained may be represented as a transfer matrix relating 
three components of the hub shear loads versus three modes of blade actuation signal 
as shown in Fig. 6-8. Note that the vertical component versus the collective mode 
actuation (the uppermost and leftmost one) is only drawn in different scale from the 
others for an easy examination of the plots. 

This transfer matrix leads to the idea of a multiple-component closed-loop con- 
troller which combines the three modes of blade actuation simultaneously. It may be 
designed for elimination of either single component of hub shear force or multiple- 
components of them at a time. For example, suppose the case of eliminating only 
the vertical components of 4P hub shear vibratory loads. From the analytical results 
presented in Fig. 5-8, the baseline amplitude for this component of the hub vibratory 
load was predicted as approximately 31 N. However, from Fig. 6-8, the amplitudes 
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Chapter 7 


Closed-loop Controller for 
Vibration Reduction in Forward 
Flight 

7.1 Overview 

The results from system identification in Chapter 6 are utilized here to design a 
closed-loop controller for hub shear vibratory load reduction. As briefly mentioned 
in Section 5.3.3., the so-called T matrix approach has been traditionally used by 
researchers investigating higher harmonic control to identify transfer matrix [8, 13]. 
This approach resulted in a rather complicated structure of a closed-loop controller in 
order to implement the modulation and demodulation phases of the scheme. However, 
it was recently found that such a complex controller structure can be reduced to a 
classical disturbance rejection algorithm [9, 10]. This results in a simple LTI feedback 
compensator. In this thesis, a closed-loop controller based on this simple structure is 
attempted. Stability of the closed-loop system is checked first. Then, the vibratory 
load reduction capability of the designed controller is demonstrated numerically by 
combining with the time domain analysis built for the ATR system. 
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d(t) 



Figure 7-1: Block diagram of the higher harmonic control system adopted by re- 
searchers using quasisteady helicopter plant model [9] 


output LTI system for a given plant transfer function G(s) [9, 10]. According to 
this derivation, the block diagram shown in Fig. 7-1 is equivalent to an LTI feedback 
compensation structure, as illustrated in Fig. 7-2. The feedback compensator, K(s), 
becomes 

(?3) 

A = Real {c^)}' B = - Imag {5UMi)} 

where T is the blade passage period. Furthermore, this LTI feedback compensator 
structure turns out to be essentially the same as a classical disturbance rejection 
algorithm, which is to eliminate an almost pure harmonic signal at constant frequency 
uj 0 with the following compensator 


H{s) = 


£i£ + C 2 

S 2 +W 2 


(7.4) 


Stability and performance issues of the closed-loop system associated with this feed- 
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Bode plot of K(s)G 0 (s) 




Figure 7-3: Bode plot of the loop transfer function 

an infinite weighting at 4P frequency, the gain of the loop transfer function goes to 
infinity at the same 4P frequency, which makes it exceed unity over a narrow region 
centered at this frequency. Then, the variation of the phase within this region will 

determine the stability of the closed-loop system. 

Another useful method of checking the stability of the controller is to examine 
the magnitude versus phase plot (Nichols plot) of the loop transfer function. Nichols 
plot of the present closed-loop system without any modification applied to K(s) is 
displayed in Fig. 7-4 (a). The stability of the system is ensured if no encirclements of 
the critical point (unity magnitude at 180° of phase). This critical point is designated 
by a small circle in Fig. 7-4. In the same plot, contours of constant disturbance 
attenuation (or amplification) are also plotted according to the following relation 

y = I (7.5) 

d 1 + K(s)G 0 (s) 

The closed contours around the critical point with positive figures represent degrees of 
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vibration amplification. The inverted U-shaped contour with thicker line represents 
0 dB boundary, where no vibration attenuation or amplification is obtained. The 
other contours indicate how much attenuation results (in dB) for the corresponding 
loop gain. Also, along the loop gain line, corresponding frequencies are designated by 
asterisks to see which frequency content is involved with a possible instability near 
the critical point. According to Fig. 7-4 (a), the present controller with the original 
K(s) turns out to have gain margin of approximately 5.4 dB, and phase margin of 
36°. In principle, this level of stability margin is regarded as enough for a general 
feedback compensator. 

7.3.2 Modified Feedback Controller 

When the gain of the controller is increased, i.e, the closed-loop gain line is shifted up- 
ward, there appears a chance of instability since it makes the gain margin diminished. 
Instability induced from this lack of gain margin is practically manifested during the 
numerical simulation, which will be described in the following section. Therefore, a 
modification on the original controller is needed to avoid this type of instability. The 
solution for tins is to alter the closed-loop gain to have new phase characteristics, 
which is shifted by 45° from its original one. This modification is implemented by 
changing the transfer function for K(s) in the following way 

K( \ _ 2fe(A / s + BW) (7.6) 

s 2 + (btt) 2 

where 

A! = Real {{A - j B)^- 45 ^} 

B' = — Imag |(A — jf?)e^ -45 

where A and B are the same parameters as defined in Eq. (7.3). The resulting 
controller with the modification generates a new Nichols plot which is shown in Fig. 7- 
4 (b). Examining the result, the phase of the modified closed-loop gam is shifted 
as desired. Notice, however, that the magnitude of the closed-loop gain is slightly 
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all the components of the harmonic transfer functions. In this section, a numerical 
demonstration of this situation is attempted using the time domain analysis. 

The resulting LTI feedback compensator, whose transfer function in continuous- 
time is described by Eq. (7.7), can be easily incorporated in the time domain analysis 
by transforming it into a discrete-time state-space approximation. Using the c2d 
command based on a Tustin approximation, which is provided in MATLAB, an ap- 
proximately equivalent discrete time state-space representation of the same feedback 
compensator can be obtained with the specific time step size selected for the analysis. 
The equivalent LTI feedback compensator which is transformed into the discrete-time 
state-space version can be expressed as 

x k+1 = Axfc + Bu fc (7.8) 

y k = Cxfc + Du fc 

where x k is the state variables internal to the controller estimated at time step k. u k 
is a sensor measurement from the plant, i.e. , hub vertical shear load measurement in 
this case, and y k is a controlled quantity of electric field applied at the blade integral 
actuators at time k in collective mode. Since the resultant controller is described as 
a third-order system in Eq. (7.7), its state-space representation involves three state 
variables, meaning that x k becomes a 3x1 vector. 

After incorporating the matrix representation described in Eq. (7.8) into the time 
domain analysis, a simulation is executed in order to evaluate the performance of the 
closed-loop controller. It is still required to establish the steady-state equilibrium for 
a specific flight condition before engaging the vibration minimizing controller. Other- 
wise, huge hub vibratory loads which are induced during the transient period before 
the steady state equilibrium is reached may generate unrealistically large control sig- 
nal. Therefore, the initial 3-second period of no actuation is applied again to obtain 
the trim condition. Then, the controller is engaged with a different magnitude of its 
gain constant discretely adjusted at 0.5, 1.0, and 2.0 within each 2-second period. 

The result of the time domain analysis is displayed as a time history of the hub 


135 


Hub Vertical Vibratory Load* 



Figure 7-6: Power spectral desnity distribution of the hub vertical shear loads for 
/x = 0.140, as = — 1°, and CV = 0.0066 without and with the closed-loop controller 

engaged 


that the vibratory component at 4P frequency is decreased by different degree for 
each gain constant applied. According to this result, more than 90% of the 4P vibra- 
tory load existing in the original response is eliminated by the controller engagement. 
However, outside the narrow band of 4P frequency, it is observed that undesirable 
increase of vibratory load components are caused by the controller. Also, another 
significant vibration amplitudes are found at integral multiples of 4P frequency (e.g., 
amplitude at 8P frequency in Fig. 7-6). These vibration components are rarely af- 
fected by the current controller design, and explain the remnant vibration after the 
4P component is eliminated. Therefore, it is recommended to improve the present 
controller into one that can eliminate disturbance in multiple harmonic components 
[25] in order to reduce other vibration components simultaneously. 
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Chapter 8 


Conclusions and Recommendations 


8.1 Summary 

T his thesis presents the numerical and experimental investigations of integral twist 
actuation of helicopter rotor blades as a mean to reduce vibration induced in for- 
ward flight. Active rotor blade with embedded anisotropic piezoelectric actuators 
is designed using the analytical framework established in this thesis. An analytical 
framework is proposed to investigate integrally-twisted helicopter blades and their 
aeroelastic behavior during different flight conditions. This framework is composed 
of frequency domain analysis for hover, and time domain analysis for forward flight. 
Both analyses are based on the same three-dimensional electroelastic beam formu- 
lation with geometrical-exactness, and are coupled with the appropriate finite-state 
dynamic inflow aerodynamics model. A prototype blade is manufactured and bench 
tests are conducted to confirm its structural characteristics and adequacy of the im- 
plementation process. Hover test of a four-bladed fully-articulated rotor is performed 
using the prototype blade with the other dummy blades. A good correlation is ob- 
tained regarding the control sensitivity of the active blade at hover condition. Based 
on the prototype blade design, four ATR test blades are manufactured in order to 
compose a fully-active-blade rotor system. This rotor system is tested in the NASA 
Langley Transonic Dynamics Tunnel in a first-ever forward flight condition for an 
integrally- twisted rotor. 
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and on the classical disturbance rejection algorithm. A preliminary stability check was 
conducted, and the resultant LTI feedback compensator was demonstrated numeri- 
cally by incorporating it in the time domain analysis model. Even with the original 
plant with all the LTP components, the designed controller showed a satisfactory 
vibration reducing capability. 

This thesis has addressed several issues related with helicopter vibration reduction 
technology with integrally-twisted rotor blades, which are: 

1 . Development of an analysis framework for the rotor system with actively-twisted 
blades in hover and forward flight 

2. Investigation of dynamic characteristics of ATR blade due to its integral actu- 
ation at bench and in hover 

3. Investigation of impact upon the ATR fixed- and rotating-system loads due to 
blade twist actuation in forward flight 

4. Estimation of harmonic transfer functions of the ATR system in forward flight 
based on the linear time-periodic system theory 

5. Design and numerical demonstration of closed-loop controller which minimizes 
helicopter vibration by integral blade twist actuation 

8.2 Conclusions 

A number of conclusions can be drawn from the current helicopter vibration reduction 
study, and they are summarized as follows. 

• Experimental structural characteristics of the prototype blade compare well 
with design goals, and predictions by the established framework correlate well 
with bench test results. 

• The design strategy of the present ATR blade prove to be appropriate by con- 
firming the actuation capability and structural integrity exhibited by the pro- 
totype blade and the modified test blades. 
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and the experimental data acquired in the amplitude of the fixed- and rotating- 
system loads. Since the discrepancy is considered to result from the less sophis- 
ticated aerodynamics model adopted, analysis with more sophisticated model 
is recommended. 

• The discrepancy in the in the amplitude of the fixed- and rotating-system loads 
also results from incomplete representation of the structural components ex- 
isting in the testing apparatus. Thus, further dynamic characterization of the 
ARES testbed is recommended. 

• System identification and closed-loop controller design was conducted based on 
the established analytical model due to the limited availability of the experi- 
mental data. However, in order to apply a closed-loop controller in the practical 
system, experimental data is recommended for final tuning of the controller. 

• Further improvement of the vibration reduction capability through optimization 
study of the active blade design parameters is suggested. This includes tailoring 
of active and passive plies included in blade lay-up, sizing of the active region, 
etc. Basic dynamic characteristics of the blade such as the natural frequencies 
may be another influencing factor for control authority in specific frequencies. 

• It is desirable to use the concept of integral blade twist actuation as a means 
of acoustic noise reduction for the rotor system. In order to analyze the active 
rotor system for such purpose, more refined aerodynamics model is required to 
accurately predict the pressure distribution around the blades in a specific flight 
condition. For a wind-tunnel experiment, quiet environments are required so 
that the reliable acoustic data may be collected. A noise-reducing closed-loop 
controller can be designed and tested by conducting the procedures suggested 
in this thesis. 
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Appendix A 

State-space Formulation for Hover 
Analysis 


Explicit expressions of the state space representation for hover analysis, Eq. (2.32), 
are provided. Since detailed expressions of the sub-matrices included in the matrices 
E and A in Eq. (2.33) are found in the appendices of [63], non-zero sub-matrices in 
the matrix B are only described in this section. 
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Appendix B 


Time Integration Formulation for 
Forward Flight Analysis 


The energy decaying time integration of the beam formulation mentioned in Sec- 
tion 2.5.2 is summarized as follows. The displacement-based form of beam formula- 
tion, Eq. (2.22), can be expressed conveniently as 

(Rp B )‘ + U[h]Rp B - (RffiX - U[< + KRfe = q ( B1 ) 


where 


U[.] = 


C B b C ba 0 

0 c Bb c ba 

0 0 
[•] 0 


Note that q in Eq. (B.l) now represents all the terms including the actuation forces 
at the right hand side of Eq. (2.21). A time discretization scheme can be applied first 
over the time step, from t~ to t~ +1 : 


B n +lPBn+l ^nPfln 

A t 


+ U 


^n+l ~ \ , 

At 


PBn + Pi?n+1 
2 


(B.2) 
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Appendix C 


AFC Distribution in the ATR 
Prototype Blade 


The following table presents the electromechanical coupling coefficients (d n and d 12 ) 
of the individual AFC packs embedded in the ATR prototype blade along its spanwise 
location. It is noted that the distribution of the packs is designed so that those with 
the best performance can be located at the inboard of the blade and the performance 
difference between top and bottom surfaces at the same spanwise location is to be 
minimized. 
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Appendix D 


Material Properties of the ATR 
Test Blade Constituents 


The following two tables present the material properties of the constituent materials 
used in the ATR test blades manufacturing. 


Table D.l: Material properties of the constituents in the ATR test blades 



E-Glass 

fabric 

EA9628 

adhesive 

Rohacell foam 
spar 

Rohacell foam 
fairing 

Thickness (/xm) 

114.3 

101.6 

- 

" 

Density (kg/m 3 ) 

1,716 

1,163 

75 

35 

E l (GPa) 

19.3 

2.38 

0.0896 

0.035 

G lt (GPa) 

4.14 

0.69 

0.0296 

0.0138 


Fiber glass 

Flexible 

Front ballast 

Strain gauge 


uni-tape 

circuit 

weight (tungsten) 

wires 

— | 

Dimension 

thickness 

width 

diameter 

diameter 


203.2 /xm 

6.604 mm 

4.7625 mm 

40x0.381 mm 

Density (kg/m 3 ) 

1,799 

3,044 

19,100 

8,900 

E l (GPa) 

48.2 

- 

- 

” 

Glt (GPa) 

5.7 

- 

- 

“ 
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Appendix E 


LTP System and its Identification 


What follows is a summary of the identification technique for LTP systems presented 
in [62]. 

E.l Characteristics of LTP system 

Linear time periodic system can be represented in a state-space form as 

x(t) = A(t)x(t) + B(t)u(t) 

y(t) = C(t)x(t) + D(t)u(t) (E.l) 

where the matrices A{t), B{t ), C(t), and D(t) are in general periodic, with period 
T. When a sinusoidal signal excites an LTP system, the system responds with the 
superposition of sinusoids not only of the input frequency uj, but also of several other 
frequencies, u) + nu p , each with its own amplitude and phase, where n is an integer, 
and u) p is the system major frequency, given as 

u> p = 27 r/T (E.2) 

The frequencies u + nuJ v are shifted harmonics, and they are often referred to simply 
as “harmonics.” 
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terms of their Fourier series, and using the harmonic balance approach, one obtains 


s n x n = ^ A n - m x m + ^ A n - m u m 

ro£ Z m£Z 

y n = ^ C n _ m X m + Dn-mUm (E.8) 

m£Z m£Z 

The s umma tions in the previous expression, Eq. (E.8), can be transformed into matrix 
form as 


sX = (A — M)X + BU 

y = CX + VU (E.9) 

The updated state vector X in Eq. (E.9) represents the original ones at various 
harmonics of the system major frequency as 




X-2 

X-l 

Xo 

Xi 

X 2 


(E.10) 


The updated system dynamics matrix A in Eq. (E.9) is a doubly-infinite Toeplitz 
matrix given by 


A = 


Aq A-i A —2 ■ ■ ■ 
Ai Aq A-i . . . 
A 2 Ai Aq 


(Ell) 


where the submatrix A n represents the n-th Fourier coefficient of A(t). B, C, and 
V can be defined in the similar manner. The frequency modulation matrix J\f is an 
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Figure E-l: LTP system block diagram with three transfer functions [62] 

Gi, and G_i, respectively. Thus, 

Y ( ju ) = Go[juj)U(juj) + G\{ju))U{jio — jujp) + G-i{ju))U{juj + ju p ) (E-14) 

In time domain, Eq. (E.14) can be expressed using the convolution as follows 

y(t) = g 0 (t) * u(t) + gi(t) * (ti(«)e**] + g-iit) * {u^e-^} (E.15) 

A linear system represented by both Eqs. (E.14) and (E.15) is depicted in the block 
diagram in Fig. E-l. 

However, since there is only one equation available in order to estimate three 
transfer functions G 0 , G u and G_i, the identification problem becomes underdeter- 
mined. This leads to the need of three different input applications for composing 
three independent equations, each of which is similar to Eq. (E.14). Due to the peri- 
odic nature of the system under consideration, it is extremely important to account 
for the time of application of each input relative to the system period T. In order 
for the system behavior to be completely analyzed, multiple identical input signals 
are applied which are evenly located over the system period. In Fig. E-2, an example 
of the input signals are shown. There, three input signals are created in sine-sweep 
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Figure E-3: Delayed input signal and corresponding output of LTP system [62] 
Then, Eq. (E.17) becomes 

1 / ^ 

Y 0 U(juj) U(juj-ju p ) U(ju + ju p ) G 0 

< Yi /3 * = U{joj) U(ju-ju p )W U{ju + jUp)W- 1 < G\ > (E.19) 

V 2/3 U{juj) U(ju> — jui p )W 2 U {juj + jWp)W- 2 - G- 1 J 


When the W terms are separated, Eq. (E.19) becomes 



Notice that the first matrix multiplied to give U in Eq. (E.20) is the one that com- 
monly arises in discrete Fourier transform. Eq. (E.20) can be simply written as 

Y - UG (E.21) 
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y and ip are similarly defined. If the transfer functions of the system as many as n h 
need to be identified, an n h xn matrix U is constructed according to Eq. (E.20) with 
an appropriately modulated and Fourier transformed vector u at each row, so that 


^{[e m ^‘ni e mj,pn u n }} 


U = 


F{[e Wl it! 

JF{[e- Wl Ui 


F{{ e ~m^ Ul 


e aj ^un}} 


e~ lj ^u n }} 


e -mjVn Un ]} 


(E.24) 


where m = ^ 


U = 


u(w 


Eq. (E.24) can be described in a more compact notation as 
— mu) p ) ... u(w — u p ) u(o>) u(w + u;p) ... u(o; + mojp) 


■\T 


(E.25) 

Similarly, Y can be constructed as the discrete Fourier transform of the vector y as 


Y — !F 


Vi 2/2 2/3 


•■■*]} 


(E.26) 


Recalling that the empirical transfer function estimate (ETFE) of a linear time in- 
variant (LTI) system involves the power and cross spectral densities of input and 
output, these spectral densities can be defined in a similar manner for the present 
LTP system as 


$uu = U* T U 

* UY = U* t Y (E.27) 
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curvature of the transfer functions, so that 


J = min [(Y - U t G) 2 + a(D 2 G) 2 ] (E.32) 

where D 2 is a second-order difference operator, and a is a weighting factor. By 
taking the derivative of J with respect to G in Eq. (E.32) and setting it to zero, the 
minimizing G can be found as 

G = [U T U + aD 4 ] _1 U r Y (E.33) 

where D 4 = D 2 D 2 . Eq. (E.33) is the final form that is utilized in the following 
system identification, and more issues on the practical implementation of Eq. (E.33) 
and solutions are provided in [62]. 
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